Supervised learning in high dimension: R
WANDIYA James, GAOLATLHE Angelah, AKAKPO Sara
Linear Regression: Pollution Example
Before beginning the analysis, it is helpful to start by visualizing the data.
| Id | maxO3 | T9 | T12 | T15 | Ne9 | Ne12 | Ne15 | Vx9 | Vx12 | Vx15 | maxO3v | wind | rain |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 20010601 | 87 | 15.6 | 18.5 | 18.4 | 4 | 4 | 8 | 0.69 | -1.71 | -0.69 | 84 | Nord | Sec |
| 20010602 | 82 | 17.0 | 18.4 | 17.7 | 5 | 5 | 7 | -4.33 | -4.00 | -3.00 | 87 | Nord | Sec |
| 20010603 | 92 | 15.3 | 17.6 | 19.5 | 2 | 5 | 4 | 2.95 | 1.88 | 0.52 | 82 | Est | Sec |
| 20010604 | 114 | 16.2 | 19.7 | 22.5 | 1 | 1 | 0 | 0.98 | 0.35 | -0.17 | 92 | Nord | Sec |
| 20010605 | 94 | 17.4 | 20.5 | 20.4 | 8 | 8 | 7 | -0.50 | -2.95 | -4.33 | 114 | Ouest | Sec |
| 20010606 | 80 | 17.7 | 19.8 | 18.3 | 6 | 6 | 7 | -5.64 | -5.00 | -6.00 | 94 | Ouest | Pluie |
Correlation
## corrplot 0.95 loaded
library("scatterplot3d")
scatterplot3d(ozone[,"T12"],ozone[,"Vx12"],ozone[,"maxO3"],type="h", pch=16,box=FALSE,xlab="Température à 12h",ylab="Vent",zlab="Ozone")
## When the points are high, this means that the rate of ozone is higher
in the air. We can notice that the ozone increases as the temperature is
increases also.
Data
## Id maxO3 T9 T12 T15 Ne9 Ne12 Ne15 Vx9 Vx12 Vx15 maxO3v
## 1 20010601 87 15.6 18.5 18.4 4 4 8 0.6946 -1.7101 -0.6946 84
## 2 20010602 82 17.0 18.4 17.7 5 5 7 -4.3301 -4.0000 -3.0000 87
## 3 20010603 92 15.3 17.6 19.5 2 5 4 2.9544 1.8794 0.5209 82
## 4 20010604 114 16.2 19.7 22.5 1 1 0 0.9848 0.3473 -0.1736 92
## 5 20010605 94 17.4 20.5 20.4 8 8 7 -0.5000 -2.9544 -4.3301 114
## 6 20010606 80 17.7 19.8 18.3 6 6 7 -5.6382 -5.0000 -6.0000 94
## wind rain
## 1 Nord Sec
## 2 Nord Sec
## 3 Est Sec
## 4 Nord Sec
## 5 Ouest Sec
## 6 Ouest PluiePCA
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Temperature and wind components contribute in 56.45% in the first
dimension, while the cloud cover contribute at 17.16% in the second
dimension.
When the wind blows in the west direction, there is a cloud cover.
Example of pollution
We consider a linear model between ozone, wind, and temperature.
modele1 <- lm(maxO3~Vx12+T12,data=ozone)
summary(modele1)##
## Call:
## lm(formula = maxO3 ~ Vx12 + T12, data = ozone)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.08 -11.99 -1.20 10.67 45.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.4242 9.3943 -1.535 0.12758
## Vx12 2.0742 0.5987 3.465 0.00076 ***
## T12 5.0202 0.4140 12.125 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.75 on 109 degrees of freedom
## Multiple R-squared: 0.6533, Adjusted R-squared: 0.6469
## F-statistic: 102.7 on 2 and 109 DF, p-value: < 2.2e-16anova(modele1)## Analysis of Variance Table
##
## Response: maxO3
## Df Sum Sq Mean Sq F value Pr(>F)
## Vx12 1 16367 16367 58.338 9.074e-12 ***
## T12 1 41244 41244 147.010 < 2.2e-16 ***
## Residuals 109 30580 281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Only the wind component and temperature at 12PM are significant. The intercept is not significant.The model only explain about 64.7% of the variance which is not good enough.
Model without intercept
##
## Call:
## lm(formula = maxO3 ~ -1 + Vx12 + T12, data = ozone)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.135 -12.914 -1.394 10.020 48.525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Vx12 2.4412 0.5523 4.42 2.32e-05 ***
## T12 4.3967 0.0811 54.22 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.85 on 110 degrees of freedom
## Multiple R-squared: 0.9688, Adjusted R-squared: 0.9682
## F-statistic: 1708 on 2 and 110 DF, p-value: < 2.2e-16When we remove the intercept which was previously not significant, the new model explain about 96% of the variance which is good.
Without wind
##
## Call:
## lm(formula = maxO3 ~ Ne12, data = ozone)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.020 -14.487 -5.115 12.571 66.470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 130.0201 4.9807 26.105 < 2e-16 ***
## Ne12 -7.9150 0.9042 -8.753 2.77e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.74 on 110 degrees of freedom
## Multiple R-squared: 0.4106, Adjusted R-squared: 0.4052
## F-statistic: 76.62 on 1 and 110 DF, p-value: 2.769e-14The adjusted R-squared explains only 40% of the variance which is very low.
##
## Call:
## lm(formula = maxO3 ~ Vx9 + Vx12 + Vx15 + Ne9 + Ne12 + Ne15 +
## T9 + T12 + T15 + maxO3v, data = ozone)
##
## Residuals:
## Min 1Q Median 3Q Max
## -53.566 -8.727 -0.403 7.599 39.458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.24442 13.47190 0.909 0.3656
## Vx9 0.94791 0.91228 1.039 0.3013
## Vx12 0.03120 1.05523 0.030 0.9765
## Vx15 0.41859 0.91568 0.457 0.6486
## Ne9 -2.18909 0.93824 -2.333 0.0216 *
## Ne12 -0.42102 1.36766 -0.308 0.7588
## Ne15 0.18373 1.00279 0.183 0.8550
## T9 -0.01901 1.12515 -0.017 0.9866
## T12 2.22115 1.43294 1.550 0.1243
## T15 0.55853 1.14464 0.488 0.6266
## maxO3v 0.35198 0.06289 5.597 1.88e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.36 on 101 degrees of freedom
## Multiple R-squared: 0.7638, Adjusted R-squared: 0.7405
## F-statistic: 32.67 on 10 and 101 DF, p-value: < 2.2e-16Only the cloud cover(Ne9) at 9AM and previous max ozone(max03v) are significant and the model explains 74% of the variance. Even though this is a good adjusted squared, the model is very complex as it has a lot of variables, most of which are not significant.
Exhaustive search
The “leaps” package allows you to perform an exhaustive search using the “regsubsets” function. For more than 50 variables, the option “really.big=TRUE” must be added to avoid excessive computation time.
## Subset selection object
## Call: regsubsets.formula(maxO3 ~ Vx9 + Vx12 + Vx15 + Ne9 + Ne12 + Ne15 +
## T9 + T12 + T15 + maxO3v, int = T, nbest = 1, nvmax = 11,
## method = "exh", data = ozone)
## 10 Variables (and intercept)
## Forced in Forced out
## Vx9 FALSE FALSE
## Vx12 FALSE FALSE
## Vx15 FALSE FALSE
## Ne9 FALSE FALSE
## Ne12 FALSE FALSE
## Ne15 FALSE FALSE
## T9 FALSE FALSE
## T12 FALSE FALSE
## T15 FALSE FALSE
## maxO3v FALSE FALSE
## 1 subsets of each size up to 10
## Selection Algorithm: exhaustive
## Vx9 Vx12 Vx15 Ne9 Ne12 Ne15 T9 T12 T15 maxO3v
## 1 ( 1 ) " " " " " " " " " " " " " " "*" " " " "
## 2 ( 1 ) " " " " " " " " " " " " " " "*" " " "*"
## 3 ( 1 ) " " " " " " "*" " " " " " " "*" " " "*"
## 4 ( 1 ) "*" " " " " "*" " " " " " " "*" " " "*"
## 5 ( 1 ) "*" " " "*" "*" " " " " " " "*" " " "*"
## 6 ( 1 ) "*" " " "*" "*" " " " " " " "*" "*" "*"
## 7 ( 1 ) "*" " " "*" "*" "*" " " " " "*" "*" "*"
## 8 ( 1 ) "*" " " "*" "*" "*" "*" " " "*" "*" "*"
## 9 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" "*"
## 10 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"The most important variable is the temperature at 12PM(T12) and the previous max ozone(maxo3v). And the cloud cover at 9AM(Ne9) in the third position.
R2 Choice
Keep Vx9, Ne9, T12, maxO3v
plot(choix,scale="adjr2")
By BIC
Keep Intercept Vx9, Ne9, T12, maxO3v
plot(choix,scale="bic")
If p>>n?
Let us consider the data on the contents of fat, sugar, flour, and water in biscuits (Osbone et al., 1984) for 72 biscuits, available in the “ppls” package. For these 72 biscuits, the compositions in fat, sugar, flour, and water are measured using a classical approach, while the spectrum is observed over all wavelengths between 1100 and 2498 nanometers, regularly spaced every 2 nanometers. We therefore have 700 observed values, or potentially explanatory variables, for each biscuit dough sample.
Typically, this study is conducted in a very high-dimensional context with \(p >> n\). On these uncooked biscuits, absorbance is measured for wavelengths between 1100 and 2498 nanometers, regularly spaced every 2 nanometers. There are thus 700 potentially explanatory variables.
Let us consider the percentage of sugars. We have \(p = 700\) variables for \(n = 40\) individuals. The classical least squares estimator is not defined.
The classical least squares estimator is not defined.
Lectrure des données
library(ppls)
data(cookie)
# extraire le taux de sucrose et les spectres
cook = data.frame(cookie$constituents$sucrose,cookie$NIR)
names(cook)= c("sucrose",paste("X",1:700,sep=""))Data
# Response
hist(cook[,"sucrose"])
## The sucrose percentage is between 8% and 22%. The distribution is
most concentrated between 10% and 20%. This shows moderate variability
in sugar content across the biscuits samples
# spectres by sucrose intensity
coul=rainbow(20)[as.integer(as.factor
(as.integer(cook[,1]-10)))]
ts.plot(t(cook[,-1]),col=coul) 
## The spectral plot shows us many colors mean that we have to have to
use multivariate method instead of univariate method.
MCO estimation not valide
cook.app=cook[1:40,]
summary(lm(sucrose ~ ., data=cook.app))##
## Call:
## lm(formula = sucrose ~ ., data = cook.app)
##
## Residuals:
## ALL 40 residuals are 0: no residual degrees of freedom!
##
## Coefficients: (661 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.17 NaN NaN NaN
## X1 13724.43 NaN NaN NaN
## X2 -31660.99 NaN NaN NaN
## X3 25153.86 NaN NaN NaN
## X4 -13071.94 NaN NaN NaN
## X5 244.71 NaN NaN NaN
## X6 24773.95 NaN NaN NaN
## X7 -26093.74 NaN NaN NaN
## X8 12801.63 NaN NaN NaN
## X9 -4709.51 NaN NaN NaN
## X10 -4891.93 NaN NaN NaN
## X11 31636.00 NaN NaN NaN
## X12 -21089.99 NaN NaN NaN
## X13 -1917.31 NaN NaN NaN
## X14 -14403.37 NaN NaN NaN
## X15 15061.66 NaN NaN NaN
## X16 22386.40 NaN NaN NaN
## X17 -54867.14 NaN NaN NaN
## X18 287.33 NaN NaN NaN
## X19 -8562.09 NaN NaN NaN
## X20 34838.19 NaN NaN NaN
## X21 17488.40 NaN NaN NaN
## X22 -23338.41 NaN NaN NaN
## X23 20792.72 NaN NaN NaN
## X24 -32664.20 NaN NaN NaN
## X25 35744.76 NaN NaN NaN
## X26 12087.43 NaN NaN NaN
## X27 -13153.13 NaN NaN NaN
## X28 -28838.43 NaN NaN NaN
## X29 22346.68 NaN NaN NaN
## X30 -17483.16 NaN NaN NaN
## X31 28184.57 NaN NaN NaN
## X32 -25085.04 NaN NaN NaN
## X33 -6570.85 NaN NaN NaN
## X34 -8131.31 NaN NaN NaN
## X35 14008.36 NaN NaN NaN
## X36 -6829.77 NaN NaN NaN
## X37 10198.76 NaN NaN NaN
## X38 -13024.01 NaN NaN NaN
## X39 14510.81 NaN NaN NaN
## X40 NA NA NA NA
## X41 NA NA NA NA
## X42 NA NA NA NA
## X43 NA NA NA NA
## X44 NA NA NA NA
## X45 NA NA NA NA
## X46 NA NA NA NA
## X47 NA NA NA NA
## X48 NA NA NA NA
## X49 NA NA NA NA
## X50 NA NA NA NA
## X51 NA NA NA NA
## X52 NA NA NA NA
## X53 NA NA NA NA
## X54 NA NA NA NA
## X55 NA NA NA NA
## X56 NA NA NA NA
## X57 NA NA NA NA
## X58 NA NA NA NA
## X59 NA NA NA NA
## X60 NA NA NA NA
## X61 NA NA NA NA
## X62 NA NA NA NA
## X63 NA NA NA NA
## X64 NA NA NA NA
## X65 NA NA NA NA
## X66 NA NA NA NA
## X67 NA NA NA NA
## X68 NA NA NA NA
## X69 NA NA NA NA
## X70 NA NA NA NA
## X71 NA NA NA NA
## X72 NA NA NA NA
## X73 NA NA NA NA
## X74 NA NA NA NA
## X75 NA NA NA NA
## X76 NA NA NA NA
## X77 NA NA NA NA
## X78 NA NA NA NA
## X79 NA NA NA NA
## X80 NA NA NA NA
## X81 NA NA NA NA
## X82 NA NA NA NA
## X83 NA NA NA NA
## X84 NA NA NA NA
## X85 NA NA NA NA
## X86 NA NA NA NA
## X87 NA NA NA NA
## X88 NA NA NA NA
## X89 NA NA NA NA
## X90 NA NA NA NA
## X91 NA NA NA NA
## X92 NA NA NA NA
## X93 NA NA NA NA
## X94 NA NA NA NA
## X95 NA NA NA NA
## X96 NA NA NA NA
## X97 NA NA NA NA
## X98 NA NA NA NA
## X99 NA NA NA NA
## X100 NA NA NA NA
## X101 NA NA NA NA
## X102 NA NA NA NA
## X103 NA NA NA NA
## X104 NA NA NA NA
## X105 NA NA NA NA
## X106 NA NA NA NA
## X107 NA NA NA NA
## X108 NA NA NA NA
## X109 NA NA NA NA
## X110 NA NA NA NA
## X111 NA NA NA NA
## X112 NA NA NA NA
## X113 NA NA NA NA
## X114 NA NA NA NA
## X115 NA NA NA NA
## X116 NA NA NA NA
## X117 NA NA NA NA
## X118 NA NA NA NA
## X119 NA NA NA NA
## X120 NA NA NA NA
## X121 NA NA NA NA
## X122 NA NA NA NA
## X123 NA NA NA NA
## X124 NA NA NA NA
## X125 NA NA NA NA
## X126 NA NA NA NA
## X127 NA NA NA NA
## X128 NA NA NA NA
## X129 NA NA NA NA
## X130 NA NA NA NA
## X131 NA NA NA NA
## X132 NA NA NA NA
## X133 NA NA NA NA
## X134 NA NA NA NA
## X135 NA NA NA NA
## X136 NA NA NA NA
## X137 NA NA NA NA
## X138 NA NA NA NA
## X139 NA NA NA NA
## X140 NA NA NA NA
## X141 NA NA NA NA
## X142 NA NA NA NA
## X143 NA NA NA NA
## X144 NA NA NA NA
## X145 NA NA NA NA
## X146 NA NA NA NA
## X147 NA NA NA NA
## X148 NA NA NA NA
## X149 NA NA NA NA
## X150 NA NA NA NA
## X151 NA NA NA NA
## X152 NA NA NA NA
## X153 NA NA NA NA
## X154 NA NA NA NA
## X155 NA NA NA NA
## X156 NA NA NA NA
## X157 NA NA NA NA
## X158 NA NA NA NA
## X159 NA NA NA NA
## X160 NA NA NA NA
## X161 NA NA NA NA
## X162 NA NA NA NA
## X163 NA NA NA NA
## X164 NA NA NA NA
## X165 NA NA NA NA
## X166 NA NA NA NA
## X167 NA NA NA NA
## X168 NA NA NA NA
## X169 NA NA NA NA
## X170 NA NA NA NA
## X171 NA NA NA NA
## X172 NA NA NA NA
## X173 NA NA NA NA
## X174 NA NA NA NA
## X175 NA NA NA NA
## X176 NA NA NA NA
## X177 NA NA NA NA
## X178 NA NA NA NA
## X179 NA NA NA NA
## X180 NA NA NA NA
## X181 NA NA NA NA
## X182 NA NA NA NA
## X183 NA NA NA NA
## X184 NA NA NA NA
## X185 NA NA NA NA
## X186 NA NA NA NA
## X187 NA NA NA NA
## X188 NA NA NA NA
## X189 NA NA NA NA
## X190 NA NA NA NA
## X191 NA NA NA NA
## X192 NA NA NA NA
## X193 NA NA NA NA
## X194 NA NA NA NA
## X195 NA NA NA NA
## X196 NA NA NA NA
## X197 NA NA NA NA
## X198 NA NA NA NA
## X199 NA NA NA NA
## X200 NA NA NA NA
## X201 NA NA NA NA
## X202 NA NA NA NA
## X203 NA NA NA NA
## X204 NA NA NA NA
## X205 NA NA NA NA
## X206 NA NA NA NA
## X207 NA NA NA NA
## X208 NA NA NA NA
## X209 NA NA NA NA
## X210 NA NA NA NA
## X211 NA NA NA NA
## X212 NA NA NA NA
## X213 NA NA NA NA
## X214 NA NA NA NA
## X215 NA NA NA NA
## X216 NA NA NA NA
## X217 NA NA NA NA
## X218 NA NA NA NA
## X219 NA NA NA NA
## X220 NA NA NA NA
## X221 NA NA NA NA
## X222 NA NA NA NA
## X223 NA NA NA NA
## X224 NA NA NA NA
## X225 NA NA NA NA
## X226 NA NA NA NA
## X227 NA NA NA NA
## X228 NA NA NA NA
## X229 NA NA NA NA
## X230 NA NA NA NA
## X231 NA NA NA NA
## X232 NA NA NA NA
## X233 NA NA NA NA
## X234 NA NA NA NA
## X235 NA NA NA NA
## X236 NA NA NA NA
## X237 NA NA NA NA
## X238 NA NA NA NA
## X239 NA NA NA NA
## X240 NA NA NA NA
## X241 NA NA NA NA
## X242 NA NA NA NA
## X243 NA NA NA NA
## X244 NA NA NA NA
## X245 NA NA NA NA
## X246 NA NA NA NA
## X247 NA NA NA NA
## X248 NA NA NA NA
## X249 NA NA NA NA
## X250 NA NA NA NA
## X251 NA NA NA NA
## X252 NA NA NA NA
## X253 NA NA NA NA
## X254 NA NA NA NA
## X255 NA NA NA NA
## X256 NA NA NA NA
## X257 NA NA NA NA
## X258 NA NA NA NA
## X259 NA NA NA NA
## X260 NA NA NA NA
## X261 NA NA NA NA
## X262 NA NA NA NA
## X263 NA NA NA NA
## X264 NA NA NA NA
## X265 NA NA NA NA
## X266 NA NA NA NA
## X267 NA NA NA NA
## X268 NA NA NA NA
## X269 NA NA NA NA
## X270 NA NA NA NA
## X271 NA NA NA NA
## X272 NA NA NA NA
## X273 NA NA NA NA
## X274 NA NA NA NA
## X275 NA NA NA NA
## X276 NA NA NA NA
## X277 NA NA NA NA
## X278 NA NA NA NA
## X279 NA NA NA NA
## X280 NA NA NA NA
## X281 NA NA NA NA
## X282 NA NA NA NA
## X283 NA NA NA NA
## X284 NA NA NA NA
## X285 NA NA NA NA
## X286 NA NA NA NA
## X287 NA NA NA NA
## X288 NA NA NA NA
## X289 NA NA NA NA
## X290 NA NA NA NA
## X291 NA NA NA NA
## X292 NA NA NA NA
## X293 NA NA NA NA
## X294 NA NA NA NA
## X295 NA NA NA NA
## X296 NA NA NA NA
## X297 NA NA NA NA
## X298 NA NA NA NA
## X299 NA NA NA NA
## X300 NA NA NA NA
## X301 NA NA NA NA
## X302 NA NA NA NA
## X303 NA NA NA NA
## X304 NA NA NA NA
## X305 NA NA NA NA
## X306 NA NA NA NA
## X307 NA NA NA NA
## X308 NA NA NA NA
## X309 NA NA NA NA
## X310 NA NA NA NA
## X311 NA NA NA NA
## X312 NA NA NA NA
## X313 NA NA NA NA
## X314 NA NA NA NA
## X315 NA NA NA NA
## X316 NA NA NA NA
## X317 NA NA NA NA
## X318 NA NA NA NA
## X319 NA NA NA NA
## X320 NA NA NA NA
## X321 NA NA NA NA
## X322 NA NA NA NA
## X323 NA NA NA NA
## X324 NA NA NA NA
## X325 NA NA NA NA
## X326 NA NA NA NA
## X327 NA NA NA NA
## X328 NA NA NA NA
## X329 NA NA NA NA
## X330 NA NA NA NA
## X331 NA NA NA NA
## X332 NA NA NA NA
## X333 NA NA NA NA
## X334 NA NA NA NA
## X335 NA NA NA NA
## X336 NA NA NA NA
## X337 NA NA NA NA
## X338 NA NA NA NA
## X339 NA NA NA NA
## X340 NA NA NA NA
## X341 NA NA NA NA
## X342 NA NA NA NA
## X343 NA NA NA NA
## X344 NA NA NA NA
## X345 NA NA NA NA
## X346 NA NA NA NA
## X347 NA NA NA NA
## X348 NA NA NA NA
## X349 NA NA NA NA
## X350 NA NA NA NA
## X351 NA NA NA NA
## X352 NA NA NA NA
## X353 NA NA NA NA
## X354 NA NA NA NA
## X355 NA NA NA NA
## X356 NA NA NA NA
## X357 NA NA NA NA
## X358 NA NA NA NA
## X359 NA NA NA NA
## X360 NA NA NA NA
## X361 NA NA NA NA
## X362 NA NA NA NA
## X363 NA NA NA NA
## X364 NA NA NA NA
## X365 NA NA NA NA
## X366 NA NA NA NA
## X367 NA NA NA NA
## X368 NA NA NA NA
## X369 NA NA NA NA
## X370 NA NA NA NA
## X371 NA NA NA NA
## X372 NA NA NA NA
## X373 NA NA NA NA
## X374 NA NA NA NA
## X375 NA NA NA NA
## X376 NA NA NA NA
## X377 NA NA NA NA
## X378 NA NA NA NA
## X379 NA NA NA NA
## X380 NA NA NA NA
## X381 NA NA NA NA
## X382 NA NA NA NA
## X383 NA NA NA NA
## X384 NA NA NA NA
## X385 NA NA NA NA
## X386 NA NA NA NA
## X387 NA NA NA NA
## X388 NA NA NA NA
## X389 NA NA NA NA
## X390 NA NA NA NA
## X391 NA NA NA NA
## X392 NA NA NA NA
## X393 NA NA NA NA
## X394 NA NA NA NA
## X395 NA NA NA NA
## X396 NA NA NA NA
## X397 NA NA NA NA
## X398 NA NA NA NA
## X399 NA NA NA NA
## X400 NA NA NA NA
## X401 NA NA NA NA
## X402 NA NA NA NA
## X403 NA NA NA NA
## X404 NA NA NA NA
## X405 NA NA NA NA
## X406 NA NA NA NA
## X407 NA NA NA NA
## X408 NA NA NA NA
## X409 NA NA NA NA
## X410 NA NA NA NA
## X411 NA NA NA NA
## X412 NA NA NA NA
## X413 NA NA NA NA
## X414 NA NA NA NA
## X415 NA NA NA NA
## X416 NA NA NA NA
## X417 NA NA NA NA
## X418 NA NA NA NA
## X419 NA NA NA NA
## X420 NA NA NA NA
## X421 NA NA NA NA
## X422 NA NA NA NA
## X423 NA NA NA NA
## X424 NA NA NA NA
## X425 NA NA NA NA
## X426 NA NA NA NA
## X427 NA NA NA NA
## X428 NA NA NA NA
## X429 NA NA NA NA
## X430 NA NA NA NA
## X431 NA NA NA NA
## X432 NA NA NA NA
## X433 NA NA NA NA
## X434 NA NA NA NA
## X435 NA NA NA NA
## X436 NA NA NA NA
## X437 NA NA NA NA
## X438 NA NA NA NA
## X439 NA NA NA NA
## X440 NA NA NA NA
## X441 NA NA NA NA
## X442 NA NA NA NA
## X443 NA NA NA NA
## X444 NA NA NA NA
## X445 NA NA NA NA
## X446 NA NA NA NA
## X447 NA NA NA NA
## X448 NA NA NA NA
## X449 NA NA NA NA
## X450 NA NA NA NA
## X451 NA NA NA NA
## X452 NA NA NA NA
## X453 NA NA NA NA
## X454 NA NA NA NA
## X455 NA NA NA NA
## X456 NA NA NA NA
## X457 NA NA NA NA
## X458 NA NA NA NA
## X459 NA NA NA NA
## X460 NA NA NA NA
## X461 NA NA NA NA
## X462 NA NA NA NA
## X463 NA NA NA NA
## X464 NA NA NA NA
## X465 NA NA NA NA
## X466 NA NA NA NA
## X467 NA NA NA NA
## X468 NA NA NA NA
## X469 NA NA NA NA
## X470 NA NA NA NA
## X471 NA NA NA NA
## X472 NA NA NA NA
## X473 NA NA NA NA
## X474 NA NA NA NA
## X475 NA NA NA NA
## X476 NA NA NA NA
## X477 NA NA NA NA
## X478 NA NA NA NA
## X479 NA NA NA NA
## X480 NA NA NA NA
## X481 NA NA NA NA
## X482 NA NA NA NA
## X483 NA NA NA NA
## X484 NA NA NA NA
## X485 NA NA NA NA
## X486 NA NA NA NA
## X487 NA NA NA NA
## X488 NA NA NA NA
## X489 NA NA NA NA
## X490 NA NA NA NA
## X491 NA NA NA NA
## X492 NA NA NA NA
## X493 NA NA NA NA
## X494 NA NA NA NA
## X495 NA NA NA NA
## X496 NA NA NA NA
## X497 NA NA NA NA
## X498 NA NA NA NA
## X499 NA NA NA NA
## X500 NA NA NA NA
## X501 NA NA NA NA
## X502 NA NA NA NA
## X503 NA NA NA NA
## X504 NA NA NA NA
## X505 NA NA NA NA
## X506 NA NA NA NA
## X507 NA NA NA NA
## X508 NA NA NA NA
## X509 NA NA NA NA
## X510 NA NA NA NA
## X511 NA NA NA NA
## X512 NA NA NA NA
## X513 NA NA NA NA
## X514 NA NA NA NA
## X515 NA NA NA NA
## X516 NA NA NA NA
## X517 NA NA NA NA
## X518 NA NA NA NA
## X519 NA NA NA NA
## X520 NA NA NA NA
## X521 NA NA NA NA
## X522 NA NA NA NA
## X523 NA NA NA NA
## X524 NA NA NA NA
## X525 NA NA NA NA
## X526 NA NA NA NA
## X527 NA NA NA NA
## X528 NA NA NA NA
## X529 NA NA NA NA
## X530 NA NA NA NA
## X531 NA NA NA NA
## X532 NA NA NA NA
## X533 NA NA NA NA
## X534 NA NA NA NA
## X535 NA NA NA NA
## X536 NA NA NA NA
## X537 NA NA NA NA
## X538 NA NA NA NA
## X539 NA NA NA NA
## X540 NA NA NA NA
## X541 NA NA NA NA
## X542 NA NA NA NA
## X543 NA NA NA NA
## X544 NA NA NA NA
## X545 NA NA NA NA
## X546 NA NA NA NA
## X547 NA NA NA NA
## X548 NA NA NA NA
## X549 NA NA NA NA
## X550 NA NA NA NA
## X551 NA NA NA NA
## X552 NA NA NA NA
## X553 NA NA NA NA
## X554 NA NA NA NA
## X555 NA NA NA NA
## X556 NA NA NA NA
## X557 NA NA NA NA
## X558 NA NA NA NA
## X559 NA NA NA NA
## X560 NA NA NA NA
## X561 NA NA NA NA
## X562 NA NA NA NA
## X563 NA NA NA NA
## X564 NA NA NA NA
## X565 NA NA NA NA
## X566 NA NA NA NA
## X567 NA NA NA NA
## X568 NA NA NA NA
## X569 NA NA NA NA
## X570 NA NA NA NA
## X571 NA NA NA NA
## X572 NA NA NA NA
## X573 NA NA NA NA
## X574 NA NA NA NA
## X575 NA NA NA NA
## X576 NA NA NA NA
## X577 NA NA NA NA
## X578 NA NA NA NA
## X579 NA NA NA NA
## X580 NA NA NA NA
## X581 NA NA NA NA
## X582 NA NA NA NA
## X583 NA NA NA NA
## X584 NA NA NA NA
## X585 NA NA NA NA
## X586 NA NA NA NA
## X587 NA NA NA NA
## X588 NA NA NA NA
## X589 NA NA NA NA
## X590 NA NA NA NA
## X591 NA NA NA NA
## X592 NA NA NA NA
## X593 NA NA NA NA
## X594 NA NA NA NA
## X595 NA NA NA NA
## X596 NA NA NA NA
## X597 NA NA NA NA
## X598 NA NA NA NA
## X599 NA NA NA NA
## X600 NA NA NA NA
## X601 NA NA NA NA
## X602 NA NA NA NA
## X603 NA NA NA NA
## X604 NA NA NA NA
## X605 NA NA NA NA
## X606 NA NA NA NA
## X607 NA NA NA NA
## X608 NA NA NA NA
## X609 NA NA NA NA
## X610 NA NA NA NA
## X611 NA NA NA NA
## X612 NA NA NA NA
## X613 NA NA NA NA
## X614 NA NA NA NA
## X615 NA NA NA NA
## X616 NA NA NA NA
## X617 NA NA NA NA
## X618 NA NA NA NA
## X619 NA NA NA NA
## X620 NA NA NA NA
## X621 NA NA NA NA
## X622 NA NA NA NA
## X623 NA NA NA NA
## X624 NA NA NA NA
## X625 NA NA NA NA
## X626 NA NA NA NA
## X627 NA NA NA NA
## X628 NA NA NA NA
## X629 NA NA NA NA
## X630 NA NA NA NA
## X631 NA NA NA NA
## X632 NA NA NA NA
## X633 NA NA NA NA
## X634 NA NA NA NA
## X635 NA NA NA NA
## X636 NA NA NA NA
## X637 NA NA NA NA
## X638 NA NA NA NA
## X639 NA NA NA NA
## X640 NA NA NA NA
## X641 NA NA NA NA
## X642 NA NA NA NA
## X643 NA NA NA NA
## X644 NA NA NA NA
## X645 NA NA NA NA
## X646 NA NA NA NA
## X647 NA NA NA NA
## X648 NA NA NA NA
## X649 NA NA NA NA
## X650 NA NA NA NA
## X651 NA NA NA NA
## X652 NA NA NA NA
## X653 NA NA NA NA
## X654 NA NA NA NA
## X655 NA NA NA NA
## X656 NA NA NA NA
## X657 NA NA NA NA
## X658 NA NA NA NA
## X659 NA NA NA NA
## X660 NA NA NA NA
## X661 NA NA NA NA
## X662 NA NA NA NA
## X663 NA NA NA NA
## X664 NA NA NA NA
## X665 NA NA NA NA
## X666 NA NA NA NA
## X667 NA NA NA NA
## X668 NA NA NA NA
## X669 NA NA NA NA
## X670 NA NA NA NA
## X671 NA NA NA NA
## X672 NA NA NA NA
## X673 NA NA NA NA
## X674 NA NA NA NA
## X675 NA NA NA NA
## X676 NA NA NA NA
## X677 NA NA NA NA
## X678 NA NA NA NA
## X679 NA NA NA NA
## X680 NA NA NA NA
## X681 NA NA NA NA
## X682 NA NA NA NA
## X683 NA NA NA NA
## X684 NA NA NA NA
## X685 NA NA NA NA
## X686 NA NA NA NA
## X687 NA NA NA NA
## X688 NA NA NA NA
## X689 NA NA NA NA
## X690 NA NA NA NA
## X691 NA NA NA NA
## X692 NA NA NA NA
## X693 NA NA NA NA
## X694 NA NA NA NA
## X695 NA NA NA NA
## X696 NA NA NA NA
## X697 NA NA NA NA
## X698 NA NA NA NA
## X699 NA NA NA NA
## X700 NA NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 39 and 0 DF, p-value: NAA linear model doesn’t give us any significant results that we can use to say something about the predictors because it shows null p-value, null adjusted R-squared and only calculate the estimate for 39 variables + the intercept.
Extract the training and test samples
cook.app=cook[1:40,]
cook.test=cook[41:72,]Ridge regression
Ici on a les paramètres \(\beta\) en fonction de la régularisation
library(MASS)
plot(lm.ridge(sucrose ~ ., data=cook.app,
lambda = seq(0,1,0.01)))
## At lambda=0, the coefficients are large and unstable. As lambda
increases, the coefficients converge toward zero
Régression ridge avec plusieurs paramètres de régularisation
plot(lm.ridge(sucrose ~ ., data=cook.app,
lambda = seq(0,0.2,0.001)))
cook.ridge=lm.ridge(sucrose ~ ., data=cook.app,
lambda = seq(0,0.2,0.001))Choix du paramètre de régularisation par validation croisée
The \(n=40\) observations in the training sample are randomly divided into 4 subsamples of size 10.
library(pls)##
## Attaching package: 'pls'## The following object is masked from 'package:corrplot':
##
## corrplot## The following object is masked from 'package:stats':
##
## loadingsset.seed(87)
cvseg=cvsegments(nrow(cook.app),k=4,type="random")Selection of the regularization parameter by cross-validation
For \(\tau\) (here \(\lambda\)) between \(0\) and \(1\), we estimate \(\beta\) and the corresponding PRESS on each subsample. We sum the 4 PRESS values to obtain the cross-validation PRESS and choose the \(\tau\) that minimizes this PRESS.
choix.kappa=function(kappamax,cvseg,nbe=100){
press=rep(0,nbe)
for (i in 1:length(cvseg)){
valid=cook.app[unlist(cvseg[i]),]
modele=lm.ridge(sucrose~.,
data=cook.app[unlist(cvseg[-i]),],
lambda=seq(0,kappamax,length=nbe))
coeff=coef(modele)
prediction=matrix(coeff[,1],nrow(coeff),
nrow(valid))+coeff[,-1]%*%
t(data.matrix(valid[,-1]))
press=press+rowSums((matrix(valid[,1],
nrow(coeff),nrow(valid),byrow=T)-
prediction)^2)
}
kappaet=seq(0,kappamax,length=
nbe)[which.min(press)]
return(list(kappaet=kappaet,press=press))
}Creation of a function to plot residuals versus predicted Y values
plot.res=function(x,y,titre="title")
{
plot(x,y,col="blue",ylab="Residuals",
xlab="Predicted values",main=titre)
abline(h=0,col="green")
}Display of results with the optimal PRESS
# execution
res=choix.kappa(0.5,cvseg,nbe=1000)
plot(seq(0,0.5,length=1000),res$press,ylab="PRESS",
xlab="lambda") # optimal value
## Base on the PRESS plot,we discover predictive performance occurs when
the lambda is very close to 0. This means that the ridge penalty does
not improve our prediction much in this case. The model works best with
little or no regularization.
Ridge regression with the optimal parameter
kappaet=res$kappaet
cook.ridgeo=lm.ridge(sucrose~.,data=cook.app,
lambda=kappaet)
coeff=coef(cook.ridgeo)
# Computation of predicted responses
fit.rid=rep(coeff[1],nrow(cook.app))+
as.vector(coeff[-1]%*%t(data.matrix(cook.app[,-1])))Plot of predictions versus true values
# True responses versus predicted responses
plot(fit.rid,cook.app[,"sucrose"],ylab="True responses",
xlab="Predicted responses")
# The line y = x
lines(cook.app[,"sucrose"],cook.app[,"sucrose"])
## By looking at the plot, the model predicts sucrose almost prefectly
because the predicted values closely match the true responses values.
This shows an excellent in-sample fit and that the chosen ridge penalty
is appropriate.
Residuals
# Computation of residuals
res.rid=fit.rid-cook.app[,"sucrose"]
# Residuals versus predicted responses
plot.res(fit.rid,res.rid,titre="")
## The points are scattered randomly around the zero line. This means
that the model’s predictions are generally accurate and unbiaised.
Prediction on the test sample and computation of prediction error
Since the data are centered and scaled, the appropriate transformation must be applied to obtain the estimated response.
# prediction on the test sample
ychap=rep(coeff[1],nrow(cook.test))+
as.vector(coeff[-1]%*%t(data.matrix(cook.test[,-1])))Prediction on the test sample and computation of prediction error
plot(ychap,cook.test[,1],ylab="True responses",
xlab="Predicted responses")
## Prediction error
mean((cook.test[,1]-ychap)^2)## [1] 9.882725The predicted and true values on the test data follow the same trend, and this show that the model generally predict well with only moderate errors.
Comparison between Ridge and OLS coefficients
# Linear model by OLS: lambda = 0
coefflm <- coef(lm.ridge(sucrose~.,data = cook.app,
lambda=0))Continuation
# Plot of coefficients from both models
matplot(t(rbind(coeff,coefflm)),type="l",col=c(2,1))
legend("topleft", col=c(2,1), lty=c(1,2),
legend=c("OLS beta", "Ridge beta"), cex=0.8)
## The plot shows that ridge regression smooths and stabilizes the
regression coefficients compared to the OLS, helping to prevent
overfitting.
LASSO: Selection of the regularization parameter
We choose 1000 values of \(\lambda\) and search, via cross-validation (on the mean squared error = MSE, using lars), for the optimal value. The plot shows the mean squared error as a function of the constraint fraction relative to \(\lambda\), \(\frac{\lambda}{|\hat{\beta}|_1}\).
# A sequence of lambda values
library(lars)## Loaded lars 1.3lambdas <- seq(from = 0, to = 1, length = 1000)Continuation
# Optimal lambda or delta by CV using lars
mse.cv <- cv.lars(data.matrix(cook.app[,-1]), cook.app[,1],K = 4,se=F,index=lambdas,use.Gram=F)
lambdas.op <-lambdas[which.min(mse.cv$cv)]From the plot, we can say the optimal lambda which minimize the cross-validation MSE is “lambda=1”.
Prediction using the optimal parameter
cook.lasso=lars(data.matrix(cook.app[,-1]),
cook.app[,1],use.Gram=F)
fit.lasso=predict(cook.lasso,
data.matrix(cook.app[,-1]),s=lambdas.op,
mode="fraction")$fitPrediction using the optimal parameter
# Plot predictions
plot(fit.lasso,cook.app[,"sucrose"], ylab="True responses",
xlab="Predicted responses")
res.lasso=fit.lasso-cook.app[,"sucrose"]Since, the points are very closed to the line y=x, we can say that the LASSO model’s predictions are very accurate.
Residuals using the optimal parameter
# Plot predictions versus residuals
plot.res(fit.lasso,res.lasso)
## The points are scattered randomly around the zero line. This means
that the model’s predictions are generally accurate and unbiaised.
Prediction error using the optimal parameter
# Mean squared prediction error
pred.lasso=predict(cook.lasso,
data.matrix(cook.test[,-1]),
s=lambdas.op,mode="fraction")$fit
mean((pred.lasso-cook.test[,"sucrose"])^2)## [1] 2.672186Regression PCR
We choose the number of components by cross-validation.
# Mean squared errors
library(pls)
cook.pcr=pcr(sucrose~.,ncomp=28,data=cook.app,
validation="CV",segments=cvseg)
msepcv.pcr=MSEP(cook.pcr,estimate= c("train","CV"))PCR regression
Mean errors as a function of the number of components.
plot(msepcv.pcr,type="l")
ncompo=which.min(msepcv.pcr$val["CV",,])-1
ncompo## 6 comps
## 6The training MSE(in red) decreases steadily as the number of components increases because adding more components always helps fit the training data better. The plot of number of components tends to approximately 6 components, which the PCR should use to get the best predictive performance.
Regression PCR
Evolution of the proportion of variance of the variables explained by each component.
# Mean squared errors
cook.pcro=pcr(sucrose~.,ncomp=ncompo,
data=cook.app)Regression PCR
# Mean squared errors
plot(explvar(cook.pcro),type="l",main="", pin=0.5)
## The plot shows us that the spectral data is highly redundant, which
is why one component captures 83% and the others are approximately
0%.
Regression PCR
We compute the prediction error on the test sample.
fit.pcr=predict(cook.pcro,ncomp=1:6)[,,6]
res.pcr=fit.pcr-cook.app[,"sucrose"]Prediction error on the test sample
plot.res(fit.pcr,res.pcr)
ychap=predict(cook.pcro,newdata=
cook.test)[,1,ncompo]
mean((cook.test[,"sucrose"]-ychap)^2)## [1] 1.443031The residual from the PCR is : “1.443031”
Regression PLS
PLS regression with up to 28 orthogonal components and residual plot with 28 components. We choose the number of components by cross-validation.
cook.pls=plsr(sucrose~.,ncomp=28,
data=cook.app,validation="CV",segments=cvseg)
print(cook.pls)## Partial least squares regression, fitted with the kernel algorithm.
## Cross-validated using 4 random segments.
## Call:
## plsr(formula = sucrose ~ ., ncomp = 28, data = cook.app, validation = "CV", segments = cvseg)Regression PLS
plot(cook.pls)
msepcv.pls=MSEP(cook.pls,estimate=
c("train","CV"))Regression PLS
We plot the evolution of the mean squared error as a function of the number of PLS components.
plot(msepcv.pls,type="l")
## The MSE tends to decrease from 0 to 4 components and flactuate at 5
components. It stabilize around 12 components.
The PLS is more precise than the PCR because it was able to minimize the MSE with only 4 components compared to 6 components chosen by PCR.
Regression PLS
We plot the evolution of the proportion of variance of the variables explained by each component.
plot(explvar(cook.pls),type="l",main="")
## The plot shows that 2 components capture around 83% of the variance
with the PLS while the PCR got only 1 component for 83%.
PLS optimal by cross-validation
## 5 comps
## 5
##Prediction error on the test sample
PLS is better than PCR.
# mean squared error
ychap=predict(cook.plso,newdata=
cook.test)[,1,ncompo]
plot(ychap,cook.test[,1])
mean((cook.test[,"sucrose"]-ychap)^2)## [1] 1.332782